References:

#Calculating metrics

#Sources: #https://www.r-bloggers.com/how-to-calculate-landscape-metrics-for-local-landscapes/ #https://r-spatialecology.github.io/landscapemetrics/reference/list_lsm.html #https://r-spatialecology.github.io/landscapemetrics/ #https://r-spatialecology.github.io/landscapemetrics/articles/articles/guide_sample_lsm.html#extract-landscape-metrics-at-sample-points #https://r-spatialecology.github.io/landscapemetrics/reference/sample_lsm.html#details

Setting up the package with the packages we plan to use

Brings in the CSVS with the buffer as a dataframe and runs a quick test to select columns from the dataframe displaying the NDSI values per buffer per site to show the values.

Bring in the shapefile for the points of interest and examine the areas

Think about RGEE and how it could be used

The parameterized model is run on the input layers using the CAPS software, written at UMass by Brad Compton and Eduard Ene. This software produces an output grid for each metric. Metrics fall into two groups: stressor metrics (such as road traffic, invasive plants, or nutrient enrichment), and resiliency metrics (similarity, connectedness, and aquatic connectedness).

Stressor metrics measure anthropogenic stressors that reduce the integrity of a site, while resiliency metrics measure the intrinsic ability of a site to maintain its ecological integrity, despite the impact of anthropogenic stressors.

Resiliency metrics, in reflecting the currentlandscape, do take into account anthropogenic stressors such as road traffic and impervious surfaces. The three resiliency metrics are based on the ecological distance among cells computed using the ecological settings variables described in Appendix D.

Scaling. Scaled metrics and IEI are scaled from 0-1; in geoTIFFs, these grids are expressed in terms of percent (scaled 0-100). Raw metrics and settings grids are scaled in original units, unique to each grid; in geoTIFFs, these grids are scaled from 0-255. The CAPS final landcover, capsland, represents landcover classes using integer classes (see Appendix H). GeoTIFF versions are already colored appropriately; legends files for QGIS, ArcView 3.3, and ArcMap are supplied for the Arc grid.

The coordinate reference system for all data is Massachusetts mainland State Plane, NAD83.

Hydrological alterations

Imperviousness imperv Measures the intensity of impervious surface in the watershed above the focal cell, based on imperviousness and the modeled “influence value” for each cell, which is the aquatic distance from the focal cell based on a timeof-flow model. Data source: landcover, streams, flow direction, watershed resistance, percent imperviousness

Integrity Metrics

Connectedness connect Measures the disruption of habitat connectivity caused by all forms of development between each focal cell and surrounding cells as well as the “resistance” of the surrounding undeveloped landscape, as well as the similarity of surroundings. A hypothetical organism in a highly connected cell can reach a large area of ecologically similar cells with minimal crossing of “hostile” cells. This metric uses a least-cost path algorithm to determine the area that can reach each focal cell, incorporating each cell’s similarity to the focal cell. Data source: landcover, ecological settings variables

Hydrological alterations

Imperviousness

  raw: http://jamba.provost.ads.umass.edu/web/CAPS2011/tiffzips/metricsraw/imperv.zip

Resiliency Metrics

Connectedness

  raw:http://jamba.provost.ads.umass.edu/web/CAPS2011/tiffzips/metricsraw/connect.zip

Vegetation

Vegetative structure

    http://jamba.provost.ads.umass.edu/web/caps2011/tiffzips/settings/structure.zip

Development & roads

Road traffic

  raw:http://jamba.provost.ads.umass.edu/web/CAPS2011/tiffzips/metricsraw/traffic.zip

##Two ways of bringing in the Tif Files

[1] One by One

[2] All at once

###[1]

#import immperviousness raster file
imperv <- raster("../external/data/Amherst_CAPS2011/imperv.tif", crs = crs(sample_points_f))

##Plot code to test tif
# plot(imperv)
# plot(sample_points_f, col = "red", add = TRUE)
# plot(ma_crop, add = TRUE)
# st_crs(imperv)

#Crop Imperv to AOI and plot  (code commented out)
imperv_cropped <- crop(imperv, y = ma_crop)
# plot(imperv_cropped)
# plot(sample_points_f, col = "red", add = TRUE)



#import connectedness raster file

connect <- raster("../external/data/Amherst_CAPS2011/connect.tif", crs = crs(sample_points_f))

##Plot code to test tif
# plot(connect)
# plot(sample_points_f, col = "red", add = TRUE)
# plot(ma_crop, add = TRUE)
# st_crs(connect)


#Crop connect to AOI and plot  (code commented out)
connect_cropped <- crop(connect, y = ma_crop)
# plot(connect_cropped)
# plot(sample_points_f, col = "red", add = TRUE)


#import structure raster file
structure <- raster("../external/data/Amherst_CAPS2011/structure.tif", crs = crs(sample_points_f))

##Plot code to test tif
# plot(structure)
# plot(sample_points_f, col = "red", add = TRUE)
# plot(AOI, add = TRUE)
# st_crs(structure)

#Crop structure to AOI and plot  (code commented out)
structure_cropped <- crop(structure, y = ma_crop)
# plot(structure_cropped)
# plot(sample_points_f, col = "red", add = TRUE)


# #Road traffic traffic Measures the intensity of road traffic (based on measured road traffic rates) in the neighborhood surrounding the focal cell, based on a logistic function of distance.
# Data source: landcover, traffic rates

#import connectedness raster file
traffic <- raster("../external/data/Amherst_CAPS2011/traffic.tif", crs = crs(sample_points_f))
##Plot code to test tif
# plot(traffic)
# plot(sample_points_f, col = "green", add = TRUE)
# plot(AOI, add = TRUE)
# st_crs(traffic)


#Crop connect to AOI and plot (code commented out)
traffic_cropped <- crop(traffic, y = ma_crop)
# plot(traffic_cropped)
# plot(sample_points_f, col = "green", add = TRUE)



###Example plot of data
# 
# par(mar = c(0, 0, 1, 0) + .1)
# plot(structure, c, axes = FALSE, nr = 4)
# plot(sample_points_f, col = "red", add=TRUE)
# 
# #Example of plotting out a single site
# par(mar = c(0, 0, 1, 0) + .1)
# plot(structure, c, axes = FALSE, nr = 4)
# sample_points_f %>% slice(1) %>% plot(col = "red", add=TRUE)



# Loading all the files together at once into a list
f <- dir("../external/data/Amherst_CAPS2011", full.names = TRUE, pattern = ".tif")
landcover_list <- list(connect_cropped, imperv_cropped, structure_cropped, traffic_cropped)
names(landcover_list) <- gsub("tif", "", basename(f))

#[2]  all together

# f <- dir("./data/Amherst_CAPS2011", full.names = TRUE, pattern = ".tif")
# landcover <- lapply(unname(f), function(x) {
#   r <- raster(x)
#   crs(r) <- crs(pts)
#   r
# })
# names(landcover) <- gsub(".tif", "", basename(f))


#Scale the images subtracting the min values of each (the zero values) and then dividing by the max value possible in the images, 255. 
lc_scaled <-lapply(landcover_list,function(x){
  (x- cellStats(x, min)) / 255
})
lc_stack <- stack(lc_scaled)

#Plot out all of the images at once
plot(lc_stack)

Calculate focal mean stat at a weight of 3000 for each variable of interest–imperviousness, traffic, and connectedness.

#Sources for understanding focal window
# http://r-sig-geo.2731867.n2.nabble.com/Raster-package-Focal-sum-in-circles-td7587683.html
# https://gis.stackexchange.com/questions/292409/why-different-results-of-mean-calculations-with-focal-in-r-and-esri-arcgis
# https://www.timassal.com/?p=2092
# https://gis.stackexchange.com/questions/287553/how-to-use-sum-of-circular-moving-window-for-each-center-cell-with-focal-in-r
# https://www.rdocumentation.org/packages/raster/versions/3.0-12/topics/focal
# https://www.rdocumentation.org/packages/raster/versions/3.0-12/topics/focalWeight



# NOTES ON FOCAL FUNCTION
#.............................................

#focal passes over each image pixel, and multiplies those weights by each pixel value in the neighborhood, and then sums those to get the mean
# It sums the values because sum is the default value of the argument “fun” in the function focal, which is why we have not even specified the argument “fun” in Block 1
# The focal function operates on the matrix, regardless of it being circular or rectangular and assigns the resulting value to the center cell. In a binary matrix values with 1 would be operated on and those with 0, ignored. A circular window is a rectangular matrix where, values occuring outside the defined radius are 0. 


# Since we are interested in the range of the recording device and because we have assigned a consistent value across the buffers for NDSI, ANT, and BIO varying per site per metric, we are going to use a focalWeight circle to replicate the 3000m buffer and pass it over the entirety of the images for imperv, connect, and structure. This will give us a raster for each and 


#2 apply circular moving window to continuous data
###
#set the focal weight, since we are using a circle, set number to the radius of the circle (in units of CRS)
#which you can conventientyly check using the check_landscape() function
# #cell resolution is 30 x 30
#...................................................................

#Now let's apply the focalweight function to get a focal matrix for each .tif file in the landcover_list 


#1500 Buffers
fw_sum_1500 <- lapply(landcover_list,function(x){
  fw <- focalWeight(x, 1500, type='circle')
  })


#2000 Buffers
fw_sum_2000 <- lapply(landcover_list,function(x){
  fw <- focalWeight(x, 2000, type='circle')
  })

#2500 Buffers

fw_sum_2500 <- lapply(landcover_list,function(x){
  fw <- focalWeight(x, 2500, type='circle')
  })



#Runs focal at a buffer weight of 1500 m radius
radius1 <- 1500
lc1500 <- lapply(1:4, function(x) {  # x <- 1
  f <- paste0("../external/data/Amherst_CAPS2011", names(lc_stack)[x], "_", radius1, ".tif")
  r <- focal(lc_stack[[x]], w = as.matrix(fw_sum_1500[[x]]), filename = f, overwrite = TRUE)
  return(r)
})
lc1500_stacked <- stack(lc1500)  # stack

# Plotting out the results
tit1 <- "1500 (m) Buffer "
par(mfrow = c(2, 2), mar = c(.5, 0, 2, 6))
lapply(1:4, function(x){
  plot(lc1500_stacked[[x]], axes = FALSE, box = FALSE, main =  toupper(c(tit1,gsub(".tif", "", basename(f[x])))))
  plot(st_geometry(sample_points_f), add = TRUE)
})

## [[1]]
## NULL
## 
## [[2]]
## NULL
## 
## [[3]]
## NULL
## 
## [[4]]
## NULL

## [[1]]
## NULL
## 
## [[2]]
## NULL
## 
## [[3]]
## NULL
## 
## [[4]]
## NULL

## [[1]]
## NULL
## 
## [[2]]
## NULL
## 
## [[3]]
## NULL
## 
## [[4]]
## NULL

Write out the raster to use with the random forest prediction

##Extract focal values to get predictors for points

Extract values from a Raster* object at the locations of other spatial data. You can use coordinates (points), lines, polygons or an Extent (rectangle) object. You can also use cell numbers to extract values. If y represents points, extract returns the values of a Raster* object for the cells in which a set of points fall. If y represents lines, the extract method returns the values of the cells of a Raster object that are touched by a line. If y represents polygons, the extract method returns the values of the cells of a Raster object that are covered by a polygon. A cell is covered if its center is inside the polygon (but see the weights option for considering partly covered cells; and argument small for getting values for small polygons). It is also possible to extract values for point locations from SpatialPolygons.

Code to create linear correlations between NDSI and IMperviousness, connectedness, and structure

Code to create linear correlations between ANT and IMperviousness, connectedness, and Structure

Code to create linear correlations between BIO and Imperviousness, connectedness, and structure

Code to create linear r squared dataframe from the rsquared summaries above